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Theoretical methods for dealing with diffusion-controlled reactions inevitably rely on some kind 
of approximation and to find the one that works on a particular problem is not always easy. In 
here the approximation used by Bogolyubov to study weakly non-ideal Bose gas, to be refereed to 
as weakly non-ideal Bose gas approximation (WBGA), is applied in the analysis of of the three 
reaction-diffusion models (i) A + A — ► 0, (ii) A + B — > and (iii) A + A, B + B, A + B — > 
(the ABBA model). The two types of WBGA are considered, the simpler WBGA-I and more 
complicated WBGA-II. All models are defined on the lattice to facilitate comparison with computer 
experiment (simulation). It is found that the WBGA describes A+B reaction well, it reproduces 
correct d/4 density decay exponent. However, it fails in the case of the A+A reaction and the 
ABBA model. (To cure deficiency of WBGA in dealing with A+A model the hybrid of WBGA 
and Kirkwood superposition approximation is suggested.) It is shown that the WBGA-I is identical 
to the dressed tree calculation suggested by Lee in J. Phys. A 27, 2633 (1994), and that the 
dressed tree calculation does not lead to the d/2 density decay exponent when applied to the A+A 
reaction, as normally believed, but it predicts d/4 decay exponent. Last, the usage of the small no 
approximation suggested by Mattis and Glasser in Rev. Mod. Phys. 70(3), 979 (1998) is questioned 
if used beyond A+B reaction-diffusion model. 



I. INTRODUCTION 



A variety of methods have been used to study diffusion- 
controlled reactions (see, e.g., refs. [1-8] for review) rang- 
ing from simplest pair- like or Smoluchowskii approach [9] 
towards more sophisticated methods such as many parti- 
cle density function formalism [4,5] and field theory [6-8]. 
In here, we focus on the field theory approach which is 
exact but, like in any other theory, one needs approx- 
imations to solve the problem. The particular way of 
making approximate calculations will be analyzed, the 
Bogolyubov's theory of weakly non- ideal Bose gas [10], 
to be referred to in the following as Weakly non-ideal 
Bose Gas Approximation (WBGA). 

The field theory is very attractive since it offers system- 
atic way to do calculation and there are well known pro- 
cedures how to represent stochastic dynamics of reaction- 
diffusion systems as field theory. Such mapping can be 
carried out for lattice- and off-lattice models, see [6] and 
[9] respectively for description. Final result is that the 
dynamics is governed by a second quantized Hamiltonian 
H(o,at), 

^(t) = -H(a,a^(t) (1) 



and system configuration at time t can be extracted from 
state vector ^f(i). All observables can be expressed as a 
field theoretic averages over products of creation and an- 
nihilation operators a* and a. 

There is no need whatsoever to resort to the use of field 
theory as all calculations could be done without it. The 
reasons for using it are mostly practical. Field theory is 
a powerful book-keeping device. First, approximations 
involved can be made more transparent and, second, due 
to this transparency it is relatively straightforward to 
borrow approximations from other problems (and related 
forms of field theory). 

However, despite its elegance, when applied to the 
reaction-diffusion systems, field theory seems to work 
with limited success. Most often the only practical pro- 
cedure of solving field theory is perturbative calculation. 
Diffusion is set up as the zeroth order problem and re- 
action is perturbation. Problem is that in the most in- 
teresting regime, when dimension of the system is low 
(bellow some critical dimension), perturbation series di- 
verges due to infra-red divergences. Any non-field theo- 
retic approach which aims at perturbative treatment will 
suffer in a similar way. 

The infra-read divergences are normally controlled in 
the framework of the renormalization group (RG) tech- 
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nique. When applying RG, the initial density is rele- 
vant parameter and grows under renormalization process. 
Growing initial densities lead to infinities in expressions 
which have to be controlled. So far, such control exists 
only in the limited amount of cases, the A+A reaction be- 
ing the simplest example. Clearly, there is a need to avoid 
perturbative treatment (together with RG) and look for 
alternative ways of calculation. In here, focus is on the 
WBGA. 

The WBGA was used to study reaction-diffusion sys- 
tems previously [8,9,11-13]. In somewhat technical 
terms, the basic idea behind WBGA is to neglect the 
products containing three or higher number of annihila- 
tion operators having non-zero wave vector in the Hamil- 
tonian. This procedure leads to closed set of equations 
for particle density and correlation functions. Solving 
these equations amounts to resummation of infinite se- 
ries of diagrams. In this way WBGA appears as a 
non-perturbative technique which can be used to con- 
trol infra-red divergences of perturbation series. In the 
following, we will distinguish between two types of ap- 
proximations to be refereed to as WBGA-I or WBGA-II 
which result in linear or nonlinear set of equations of mo- 
tion respectively. 

On the more intuitive basis, the neglect of fluctuations 
on the shorter length scale (large k vectors) has to do 
with the slowness of diffusion. For initial conditions when 
particle density is uniform there is no k 7^ component 
in average particle density and the density profile is flat. 
Since diffusion is a very slow process one expects that 
this profile is kept all the time up to a small perturba- 
tion. This picture is likely to be true when particle den- 
sity is low when compared to the reaction range, which is 
summarized in the condition that nr d <C 1 where n and 
r denote the particle density and the reaction range [9]. 
Naturally, there are other scales in the problem (which 
change in time) and this kind of thinking can not be em- 
ployed without encountering problems [4] . In reality, one 
really has to test the WBGA on the particular model to 
see whether it works. 

To analyze properties of the WBGA method, we apply 
it in the analysis of three diffusion-controlled reactions 
(i) the A+A, (ii) the A+B and (iii) the ABBA model. 
The ABBA model includes A+A, A+B and B+B reac- 
tions simultaneously which allows for competitions be- 
tween A+A and A+B reactions (see ref. [14,15] for more 
details). The A+A and A+B models have been studied 
by variety of methods (see, e.g., ref. [16] for review). In 
here, the focus is on the WBGA studies. Apart from few 
initial studies of the A+A and A+B models the WBGA 
have not been widely used. These early applications of 
WBGA on A+A and A+B reaction-diffusion systems are 
mostly done for stationary situation (and decay from it) 
and are listed in refs. [11] (A+A) and [12,13] (A+B). 
Also, the WBGA reappears in paper [6] where it is used 
to study the A+B reaction. 



The purpose of this work is twofold. First, we study 
A+A and A+B reactions with uncorrelated initial con- 
ditions with particles Poisson-distributed [17] at initial 
time t = 0. In previous work which employed WBGA 
method in refs. [11-13] it was assumed that particles are 
initially correlated, with a initial condition being pre- 
pared in a somewhat special way [18]. Second, once the 
workings of the WBGA approximation is illuminated, 
we apply the WBGA method to analyze the more com- 
plicated ABBA model. This model was studied previ- 
ously with field theoretic RG technique using e-expansion 
and numerical simulation [14,15]. The ABBA model de- 
scribes a mixture of the A+A and A+B reactions which 
occur at the same time. The A+A (A+B) component 
bias kinetics towards d/2 (d/A) decay exponent and it is 
not a priori clear which exponent will prevail. The anal- 
ysis in [14,15] shows that both A and B particles decay 
with d/2 decay exponent and different amplitudes. 

The A+A model has been studied by variety of meth- 
ods and lot is know about it's behavior. The study most 
relevant for this work is given in ref. [19]. For uncorre- 
lated initial conditions particle density n decays asymp- 
totically as n ~ A{Dt)~ d / 2 for d < 2 while there is 
logarithmic correction at d = 2 as n ~ ln(Dt)/(8irDt), 
where d denotes dimension of system, t is time and D de- 
notes diffusion constant. Amplitude A is universal num- 
ber independent on the details of the model. Also, the 
exact so lution of the model is possible at d = 1 where 
l/V^Di [20]. 

The d/2 decay exponent of the A+A reaction is ro- 
bust in a sense that even simplest pair like approach 
(e.g. Smoluchowskii) predicts this exponent. Thus, it 
is natural to test WBGA on the A+A model. It will 
be shown that, somewhat surprisingly, WBGA fails to 
describe A+A reaction: WBGA-I predicts d/A decay ex- 
ponent while WBGA-II gives mean field exponent. 

Taking slightly more complicated reaction scheme, 
such as the A+B, leads to very hard problem where cal- 
culation methods start to differ in their predictions. The 
nature of the asymptotics of A+B reaction when initial 
number of A and B particles is equal has been hotly de- 
bated for decades. For example, Smoluchowskii approach 
predicts d/2 decay exponent [4,5,9]. Field theoretic RG 
method predicts correct d/A decay for 3 < d < 4 [21], 
while for other dimensions nothing can be said in the 
RG framework. Strangely enough, seemingly the most 
sophisticated coupled-cluster technique argues for mean 
field exponent [22]. The true decay is governed by d/A 
exponent as shown by work of Bramson and Lebowitz 
which provided strict mathematical proof [23]. It will 
be shown that when particles are initially uncorrelated, 
i.e. Poisson-distributed, both the WBGA-I and WBGA- 
II predict d/A decay exponent. Strangely enough, WBGA 
fails on simpler A+A reaction while it describes properly 
more complicated A+B reaction. 

The paper is organized as follows. In section II the 
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model to be studied is specified in detail, and mapping 
to the field theory is described. The outcome of this sec- 
tion is Hamiltonian which describes the most general two 
species reaction-diffusion model. In section III equations 
of motion for density and correlation functions are de- 
rived using WBGA method. In section IV the WBGA ap- 
proach is used to describe A+A model, with application 
of WBGA-I discussed in section IV A, and application of 
WBGA-II in section IV B. Also in section IV A equiva- 
lence of dressed tree calculation with WBGA-I is shown. 
Section IV C discusses modification of Kirkwood super- 
position approximation in the spirit of WBGA. In section 
V A+B model is studied by using WBGA-I and WBGA- 
II methods. Finally, the analysis of the most complicated 
of three models studied here, the ABBA model, is pre- 
sented in section VI. Workings of WBGA-I are analyzed 
in section VIA and of WBGA-II in section VI B. The 
step by step merger of WBGA and Kirkwood superposi- 
tion approximation is discussed in VIC. The summary 
and outline of future work is given in section VII. 



II. THE MODEL AND THE MAPPING TO THE 
FIELD THEORY 



The mapping to the field theory will be carried out 
for the most general two-species reaction-diffusion model. 
The A+A, A+B and the ABBA model will be obtained as 
special cases. It will be assumed that particles can not be 
created, neither by external source, nor by birth process. 
In the most general version of the two species reaction- 
diffusion model each of the A or B particles jumps onto 
the one of the neighboring lattice sites with rate (or dif- 
fusion constant) Da or Db, irrespectively of occupation 
number of the site where particle jumps to. Apart from 
diffusing particles annihilate in pairs. Particle p sitting 
at lattice site x and particle v sitting at the site y are as- 
sumed to annihilate with rate a^ y . This is schematically 
denoted as 



p(x) + v{y) -» 



(2) 



where p,v = A,B and x, y = 1, 2, 3, . . . , V. It is assumed 
that there are L lattice sites in one direction, and for d 
dimensions total number of sites equals V = L d . Labels 
x and y denote lattice sites. 

The stochastic dynamics of the most general two- 
species model is governed by the master equation of the 
system. The master equation describes transitions be- 
tween different configurations c. The configuration of 
the system is specified by the occupancy of lattice sites 
at certain time t; 

c = (ni,mi, . . . ,n x ,m x , . . . ,n y ,m y , . . . ,n v ,m v ) (3) 

n x and m x count number of A and B particles re- 
spectively at x-th lattice site; n x ,m x — 0,1,2,... and 



x = 1,2, ... ,V. The quantity of interest is P(c; t), which 
denotes probability that at given time t system is in the 
configuration c. 

Diffusion and reaction contribute independently to 
change of P(c; t); 



d_ 

dt 



P(c;t) = Pn(c;t) + P R (c;t) 



(4) 



The influence of diffusion is described by 



P D (c; t) = ^{DaKth + l)P(/n x + 1, n e - 1/; t) 

x e(x) 

-n,P(c;t)} + D B [(m, + 1) 

xP(/m x + l,m e -l/;t)-miP(c;t)}} (5) 

where notation P(/ . . . /; t) indicates that content of the 
sites specified by / . . . / has been modified in c. The 
X^ e ( x ) denotes sum over all first-neighbor sites of the site 
x. The change of P(c; t) driven by reactions is given by 



P R (c;t) = 5> : 



aa (n x + 2)(n x + l) 



P(/n x + 2/;t) 



■ \- rb (mx + 2) (m x + 1) 
+ Z^^xx 2 P{/m x + 2/ ,t) 

X 

+ ^2^xx(n x + l){m x + l)P(/n x + l,m x + l/;t) 

X 

+ E a xy(n x + 1)K + l)P(/n x + l,n y + 1/; t) 

x>V 

+ 2 + l)fo/ + + l,m v + l/; t) 



x>y 

+ J2 a xy( n * + l )( m v + l ) p U^x + 1, m y + V; t) 

EAA n x(n x — 1) BB m x( m x ~ 1) 

^ xx 2 / j ® xx 2 

X " X 

a xx n x m x + ^ cr xy n x n v 



x>V 



+ X! a xy m xm y + Y VxyKxmy 
x>y x =iy 



P(c;t) (6) 



The form of the master equation describing reactions can 
not be made more compact since distinction has to be 
made between reaction on the same and different lattice 
sites. For example, the term describing A+A reaction at 
the same lattice site has the pre- factor n x (n x — 1) while 
for two different sites another form n x n y has to be used. 

The master equation (4) is the first order differential 
equation and one has to specify initial condition in the 
form P(c; 0) for all c. In the following it will be assumed 
that P(c; 0) is Poisson distribution with averages denoted 
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by no, a and no,s for A and B particles respectively. This 
means that P(c; 0) factorises, 

P(c; 0) = n x =o,i,2,...,vp(n x , n , A )p{m x , n , B ) (7) 

where p(n, n) denotes Poisson distribution function 

p(n,n) = e~ n n n /n\ , n = 0,1,2,... (8) 

Eq. (4) can be translated into a Schrodinger-type equa- 
tion (1) by defining state vector, 

= 2 P{ni,mi,...,n v ,m v ;t) 

ni ,mi ,...,ny ,my 

x(aj) ni (6 t i) mi -(4) nv (^) mv |0> ( 9 ) 

where aj., a x , &£, and 6^ are creation and annihilation 
operators for A and B particles respectively; 



[a x ,a\\ = 5 XtV , [b x ,b\] = 5 x , y 
Wx,b\] = [al,b y ] = 



(10) 



and [I, II] = I II — II I denotes commutator. One can 
see that Eq. (4) is equivalent to the Eq. (1) provided dy- 
namics is governed by the second-quantized Hamiltonian 
H = Hd + Hr. Term Hd describes diffusion 

Hu = J2 2^4k ~ ° e) + Dsb ^ - b ^ ( n ) 

x e (x) 

and reactions 

Hr = IJ2 °xy A ( a H - fyxOy + 

x,y 



J2 a -y ( - a - b l- 1 ) a - b y 



(12) 



Please note that expression for i?R has the more compact 
form than the expression for related part of the master 
equation (6) as all sums in (12) over x and y are unre- 
stricted (a first sign of convenient book-keeping). Also 
the fact that reactions A+A, B+B and A+B happen in- 
dependently from each other influences the form of Hr ; 
contributions describing A+A (first line of Eq. 12), B+B 
(second line) and A+B (third line) enter additively. 

Once H has been specified ty(t) can be found from (1), 



(13) 



In general, form of depends on the initial particle 
distribution and for Poisson-distributed particles equals 



= exp 



n ,A 2(4 - 1) + n o,B ^2(bl - 1) 



|0) (14) 



Since is available (at least in principle) one can cal- 
culate various observables with following recipe 



(0(n x (t))) = (l\0(4a x )\*(t)) 



(15) 



where 0(n x ) denotes any function which depends on par- 
ticle numbers in arbitrary way. The generalization of 
Eq. (15) for more complicated form for 0(n x ,m y , . . .) is 
trivial. Vector (1| is given by 



(1| = (0|exp 



^2{a x + b x ) 



(16) 



being left eigenvector of a\ and b\ with eigenvalue 1. To 
see that Eq. (15) indeed evaluates observables correctly 
one can expand function O in Taylor series in n x and 
check that for every term in Taylor series Eq. (15) works. 

Following ref. [6], since (1|4 + 1 0, it is useful to trans- 
form Eq. (15) slightly in order to obtain form where (1| 
changes into the vacuum state (0|. This makes calcula- 
tion somewhat easier since one can use property of the 
vacuum that (0|4 = 0. To do the transformation follow- 
ing identity proves useful, 



(110(4,0^)10) = <0|O(4 + l,a*)|0> 



(17) 



and by using (17) Eq. (15) can be written in the form 

(0(n x (t))) = (0|O((4 + l)a x )e- m \^) (18) 

where bar over H and "fo indicates that substitution 
4-4 + 1, + x = 0,l,2,...,V (19) 

have been made. The is given by 



n ,A 24 + n ,B 2 b l 



|0) (20) 



and H = Hd + Hr. After the shift Hd does not change 
the form, i.e. Hu — Hd, while changes into 

^ R = \ 2 a xy( a l + a l + 44)°» a » + 

x,V 

+^20 6 * + & £ + W>A + 

x,v 

+ ^< B (4 + bl + albl)a x b v (21) 

x,y 

In the following we compactify notation and use 
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(f(a x ,al,...)) = (0\f{a x ,al,...)e- Ht \^ ) (22) 

where / denotes any function of operators a x , , b x and 
b^ y for any x and y. Also notation 



(t) = e _1 "*o 



(23) 



will be useful. 



To avoid effects of boundaries periodic boundary con- 
ditions are assumed and one can introduce Fourier trans- 
forms for operators, 

a * = 7= E e ^ . 4 = -)= E e " ifex 4 ( 24 ) 



ikxiA 



b * = jv^ eikXbk > bl = 7v^ e ~ lKXbk 



and reaction rates 



ikx 



(25) 



(26) 



with p, v = A, B. For convenience, Fourier transforms 
are defined slightly differently for operators and reaction 
rates (just to reduce explicit occurrence of V in expres- 
sions later on). Also, it is useful to express Hamiltonian 
H in terms of creation and annihilation operators in k- 
space. Starting from (11) and (21), and using (24)-(26) 
gives 

H U = D A ^2 k 2 a{a k + D B £ k 2 b\b k + 0(k 4 ) (27) 



and 



H R = -j= (Jq A a\a k - q a q 

1 X - ^ 4 4 + + 

+ ^?E CT « ala\a k ^ q a l+q 
q,k,l 

+-^E< b& ^-a 

Q-k 

+-i 51 a £ B ( a l a k-qb q + b{b k - q a q ) 
(j.fe 

+17 E a i B a k b W-q h l+ q 



(28) 



q,k,l 



Please note that Eq. (27) is an approximation which is 
valid for small k. The form in Eq. (28) is exact. Once the 
explicit form of H is known one can proceed with calcu- 
lation of particle density, as shown in the next section. 



III. EQUATIONS OF MOTION FOR DENSITY 
AND THE CORRELATION FUNCTIONS 



Equation (18) summarizes how observables are calcu- 
lated within the field theory formalism. We continue with 
the specific case of local particle densities n A (x,t) and 
n B (x,t), which can be calculated with O — a x a x and 
O = b x b x respectively; using (18) one gets n A (x, t) = (a x ) 
and n B (x,t) = (b x ), or more explicitly 

n A (x, t) = (OKI*) = (0Ke- fit |* o ) (29) 
n B (x, t) = (0|6 X |*) = (0|^e- St |§o> (30) 



If initial conditions are translationally invariant, which 
is the case for Poisson-distributed particles, local parti- 
cles densities are position independent, i.e. n A (x, t) = 
n A (t) and n B (x, t) = n B (t). Also at t = n p (0) = 
n , p with p = A,B. It is convenient to re-express 
n A (t) and n B (t) as n A (t) = yJ2 x n A{x,t) and n B (t) = 
yJ2 x n B{ x ,t)- Using field theory n A (t) and n B (t) can 
be calculated from 



2 i 

n B (t) = ±(J2b x ) = -^=(b ) 



(31) 
(32) 



where the sum over x divided by \/V was recognized as 
k = component of a k and b k ■ 

The equations of motion for n A (t) and n B (t) can be 
derived as follows. n A (t) and n B (t) are proportional to 
(a ) and (b ) respectively, and it is sufficient to derive 
equations for them. Taking time derivative of (do) and 
(bo) gives, 



d_ 

dt 

d_ 

dt 



(a ) = -([a ,H}) 
(b ) = -([b ,H]) 



(33) 
(34) 



The commutator emerges when time derivative acts on 
exp(— Ht) [24]. Evaluating the commutator with H given 
in (27) and (28), and taking into account (31) and (32) 
gives 



dn A 
dt 



dn B 
dt 



[<j AA n A n A + a AB n A n B + 



fe/O 



+ °k l k ) 



(35) 



[a BB n B n B + a AB n A n B + 



+ yEK i 

fe/0 



BB-rBB 
1 k 



AB-nAB\ 



(36) 
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where 



1 k 


= (»/ta- 


-*> 


1 k 


= (bkb- 




r AB 
1 k 


= (a fe 6_ 


-*> 



(37) 
(38) 
(39) 



To derive (35) and (36) the sum over k is split into k = 
and k + 1 parts and assumption is made that k = com- 
ponents are non- fluctuating, i.e. {a a ) ~ (a )(a ) and 
likewise for (a &o) and (&oM (thermodynamic limit). 

The equations (35) and (36) involve correlation func- 
tions with p,v = A 1 B which we proceed to calculate. 
The time evolution of correlators is governed by [O, H] 
with O = CLkd-k, akb-k, bkb-k- Evaluating commuta- 
tors gives equations of motion, 



dt 



T AA = -2D A k*T AA [a A 



AA n\+ 



\ AAf-iAA 
' a q 1 k-q 



V ^ 

q7 tk 



+ 2K 



AA 



+ a AA )n A T AA + 



r^f) 



o^ AB ^ r AA ~\ 
2<7 n B i k J 



(40) 



d_ 

dt 



U = -2l) B k L k - [a k 



+ 



\ " BB r BB 
q^k 



2(4 



BB 



_i_ ~BB 
+ <*k . 



n B r B B + 



AB 



n B (r£ B + r_ fc 
-(D A + £> B )fc 2 r^ 



9 r AB 



o^ AB „ r BB 1 
2cr o nAl k 



(41) 



y q k-q 

q^k 



AA 



\o AB n A n B + 



<?k A )nA+ 



+ (a BB + ^ B )n B + a AB (n A + n B )] T AB - 



-a AB (n A T 



BB i tAA\ 

k +n B T k ) 



(42) 



Eqs. (40), (41) and (42) are approximate since correla- 
tors of type (aka q ai) with k ^ 0, q ^ 0, and I ^ are 
assumed to be small. This is exactly the content of the 
WBGA. 

Few comments about the inversion symmetry (in k- 
space) of T p k p, p = A, B are in order. By construc- 
tion, Y p k p {t) = r p _ p k (t) for p = A,B since these func- 
tions represent correlations for same operator types. It is 
not obvious whether this property holds for T AB which 
describes correlations for different operators type. Us- 
ing Eq. (42) to derive equation of motion for T k \t) = 
r fe S W - rA f (t) on e can see that rjf } (i) = is solution 

provided T k \o) = 0. The initial conditions for T p k with 
pv = A,B are given by 

r r(°) = hfiVn ,pn QtV , p, v = A, B 
where 8 X , V denotes Kronecker delta- function, 



(43) 




Also, notation 



(44) 



(45) 



will be used. Thus T k '(0) is indeed zero. It follows 
that V AB {t) = T AB (t) for every t > 0. In the following 
whenever T k B + T AB appears in the equations of mo- 
tion, we will use assumption of inversion symmetry and 
shorten expression to 2T AB . Naturally, if this inversion 
symmetry is broken at t = one has to keep full form 
Y AB + T AB in the equations of motion, but such case is 
not considered here. 

The last terms in Eqs. (40)-(42), which are product of 
density and correlator, and appear to be third order in 
density 0(n 3 ), lead to non- linear equations of motion. 
These terms come from averages of the type (aka q ai) 
where one of {k, q, 1} momenta is zero while remaining 
two are not. One example would be {a^auai) which is 
approximatively equal to (ao}(akai). (Same discussion 
would apply irrespectively which operator one uses in 
average, ak or b k .) Being third order in density, it is 
tempting to neglect these 0(n 3 ) terms, since one expects 
that leading contribution should come from terms which 
are second order in density 0(n 2 ), e.g. first and second 
terms on the right hand side of Eq. (40) and likewise for 
Eqs. (41) and (42). To test the effect of neglecting or 
keeping 0(n 3 ) terms two approximations will be stud- 
ied WBGA-I or WBGA-II with 0(n 3 ) terms taken away 
or kept in the calculation. WBGA-I approximation re- 
sults in a linear set of equations which makes the ana- 
lytical analysis possible. Also, as formulated above, the 
WBGA-II and WBGA appear to be equivalent. In the 
following, the term WBGA will imply both WBGA-I and 
WBGA-II. 

The equations of motion for the three models we wish 
to study, the A+A, A+B and ABBA are easily extracted 
from the most general form given in (35)- (36) and (40)- 
(42). To obtain equations of motion for the ABBA model 
one simply sets D A = D B = D and 



AA 

xy 



BB 
'xy 



X5. 



x,y ) 



AB 



55 x ,y 



(46) 



An obvious distinction is being made between Kronecker 
delta-function S x , y and the symbol S which denotes the 
reaction rate. Thus particles have to meet at the same 
lattice site in order to react. Further, to obtain the A+A 
model one simply sets 6 = (this decouples A+A and 
B+B reactions, i.e. particles A and B move and react 
independently of each other). To get A+B model one 
takes A = which rules out the A+A reaction. In the 
following section we continue with analysis of the A+A 
model within WBGA framework. 
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IV. WBGA APPLIED TO THE A+A REACTION 



Using (46) with 5 = in (35)-(36) and (40)-(42) gives 
equation of motion for the density 



9n , r 2 



and the correlator 



dt 



2Dk 2 Y k - A[n 2 + $] - 4Anr fc 



(47) 



(48) 



Letter A has been dropped on ua and to simplify 
notation, and likewise ti^ a is shortened to tiq. $(t) is 
implicitly defined by correlators, 



(49) 



fc^O 



Thus equations above are meant to describe the model 
where only one type of species, A, jumps on the lat- 
tice and particles have a chance to react only when 
at the same lattice site. The equations above will be 
solved in the following two subsections using WBGA-I 
and WBGA-II approaches. 



A. The WBGA-I approximation 

In the WBGA-I approximation, when term propor- 
tional to nTk is dropped, Eq. (48) can be studied an- 
alytically. Uncorrelated (Poisson-like) initial condition is 
described by Tk = SkfiVn-o and solution of (48), with 
k 7^ 0, reads 

T k (t) = -A f dt'e- 2Dk2 ^ [n{t'f + *(*')] (50) 
Jo 

Please note that the To(i) is determined from To(t) = 
Vn{t) 2 (thermodynamic limit) and not from (50). Sum- 
mation of Eq. (50) over k ^ and division by V gives 

$(t) = -A / dt'G(t - t')[n(t') 2 + *(*')] (51) 
Jo 



where 



(52) 



fc#0 



was introduced. As shown in the appendix A, for large 
lattice size when V — > oo, expression above can be ap- 
proximated as 



G(t-t') w [87rD(t-t' + V )r d/2 



where »7 = g^y. 

Equations (47), (51) and (53) completely specify n(t). 
It is not possible to solve them analytically, however, 
large time behavior of n(t) can be extracted. To do this 
we introduce ip = n 2 + <f> and rewrite Eqs. (47) and (51) 
as 

| = -V _ (54) 

ip(t) = n(t) 2 - A / dt'G{t - t')<p(t') (55) 
Jo 

which completely specify n(t). 

By using Laplace transform it is possible to transform 
equations (54) and (55) into a single equation. Laplace 
transform is defined as 



X(s) 



= f 

Jo 



dte- st X(t) 



(56) 



For X — n, cp same symbol will be used for Laplace trans- 
form as for the original function. The only exception to 
the rule are two cases. For X(t) = n(t) 2 , X(s) = n 2 {s), 
while for X(t) = G(t), X(s) = g(s). 

Taking Laplace transform of Eq. (55) one gets ip(s) = 
n 2 (s) — \g(s)ip(s), and combining it with p(s) = (sn(s) — 
n )/\ from (54) gives 



9 (s) + i 



[n - sn(s)] 



(57) 



The g(s) is the Laplace transform of G(t), 

g{s) = (8^)- d/2 e , ' s s d/2 - 1 r(l - d/2, r/a) (58) 

T((3,x) denotes incomplete Gamma function; 

poo 

V{fi,x) = / duu- 1+fi e- u (59) 

J X 

The analytic continuation of r(/3, x) is possible. For non- 
integer /3 and (3 = 0, T(/3, z) is multiple-valued function 
of z with a branch point at z = 0, and has no poles. 
Dividing by g(s) + 1/A Eq. (57) results in 



sn(s) - n Q = -A ff(s)n 2 (s) 



(60) 



(53) 



where A e fj(s) denotes Laplace transform of the effective 
reaction rate, 

^ = TTJgV) (61) 
Finally, taking inverse Laplace transform of (60) gives 
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dn 
~dt 



= - f dt'\ eS {t - t')n(t') 2 
Jo 



(62) 



Equations (54)- (55) and Eq. (62) are fully equivalent. 
Both equations have been dealt with before, however, in 
a very different context. Eq. (62) was obtained in ref. [19] 
through a diagrammatic technique and referred to as the 
dressed-tree calculation. Thus in here we have shown 
that WBGA-I is equivalent to the dressed tree calcula- 
tion. Eqs. (54)-(55) have been also obtained in the study 
of the entirely different A+B model in ref. [6]. 

Studies [19] and [6] suggest contradictory result: 
ref. [19] argues that the dressed tree calculation gives d/2 
decay exponent for particle density, while ref. [6] argues 
for the d/4 decay exponent. In ref. [19] it was incor- 
rectly concluded that dressed-tree calculation results in 
the d/2 exponent, which basically came from balancing 
wrong terms in Laplace transformed version of Eq. (62) . 
In here calculation done in ref. [19] will be repeated to 
show how to balance terms correctly. Also, calculation 
will justify approximations employed in ref. [6] more rig- 
orously as the method of calculation employs Laplace 
Transform and well know Tauberian theorems which re- 
late small s with large t behavior. In this way all the 
approximations are controlled. 

To extract asymptotic behavior for n(t) from Eq. (62) 
one assumes that at large times density decays as 



i(t) ^A(n + ty a 



(63) 



A and a denote amplitude and exponent of decay to be 
found, n is introduced as regulator for small t so that 
Laplace transform of n(t) and n 2 (t) exist; 



n(s) = Ae s »s a - 1 r(l-a,iJ,s) 
n 2 {s) = A 2 e s ^s 2a - 1 T(1 - 2a, (is) 



(64) 
(65) 



and please note that n(s) 2 ^ n 2 (s). To extract asymp- 
totics one inserts (64)- (65) and (61) into (60) and ex- 
pands in small s (to extract leading order behavior for 
large t) and matches the most dominant terms. 

The expansion of g(s) for small s is given by 



9(s) 



(8TrD)- d / 2 e sv x 



r(i - d/2)^/ 2 - 1 + ' - + o(s) 



d-2 



(66) 



for d ^ 2, 4, 6, For d = 2 one has 

g{s) = (SttD)-^^ [- lE - ln(r?s) + 0( V s)} (67) 

where is Euler constant. Please note that the behav- 
ior of g(s) for small s is qualitatively different for d < 2 
and d > 2 which has to do with recurrence of random 
walks bellow and above d = 2. For small s and d < 2 



g(s) oc s d l 2 ~ x while for d > 2 g(s) = const. At d = 2 
there is logarithmic dependence on s. The term e sri can 
be neglected if leading order behavior for small s (large 
t) is sought for. 

At the moment we focus on the d < 2 case. Inserting 
approximate formulas above for g(s) into (61) gives 



Aeff(s) 



(8irD) d / 2 
T(l-d/2) ! 



-d/2+1 



d < 2 



(68) 



Since value for a is not known one has to separate various 
cases: expansion for n(s) reads 



n(s) = A 



[r(l -a)s"- 1 +0(1)] a < 1 



a-l 



+ 0(s a - 1 ) 



a > 1 



(69) 



and likewise for n 2 (s) 



[r(l - 2a)s 2a - 1 +0(1)] 2a <1 
n 2 (s) =A 2 I r i- a . „ , i (70) 



2a-l 



2a > 1 



Inserting small s expansions (68)-(70) into (60) gives 

A[s a T{l - a) + 0{s)] -n = 
A 2 i^Dfl 2 



r(l-d/2) 



s 2a - d l 2 T{l-2a) + 0{s 1 - d l 2 ) (71) 



Also, please note that there arc two different forms to use 
for n{s) and n 2 (s) in (69) and (70) and the ones used in 
(71) were for a < 1 and 2a < 1 respectively (same choice 
was made in ref. [19]). Once a is found, one has to check 
these conditions on a for self consistency. There are two 
ways to match the terms in (71), (a) as in the ref. [19], 
and (b) in a way related to the work in ref. [6] . We begin 
with first case. 

Balancing s a term on the left hand side of (71) with 
s 2a-d/2 Qn j. ne rigi^ hand side, gives a = d/2 and 

A a = --sin(7rd)r(d)r(l - d/2) 2 {%TTD)- d/2 (72) 

Also from a < 1 and 2a < 1 one has constraint that 
d < 1. However, for d < 1 the term sin(7r<i) is positive 
which makes amplitude A a negative. Thus all physical 
conditions can not be met with this type of matching. 
In ref. [19] the condition d < 1 [coming from the fact 
that first row is used in (70)] was overlooked (if d > 2 is 
allowed amplitude A a is perfectly acceptable). 

The a = d/2 scenario can still turn out to be true. 
With this choice of a and d < 2 condition coming from 
(68) the second row in (70) have to be used. Again, car- 
rying out similar type of matching procedure would give 
negative amplitude. Finally, the a = d/2 avenue has to 
be given up. 
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At this stage one is left by the second (b) way of bal- 
ancing, i.e. matching constant no term on the left hand 
side of (71) with s 2a - d / 2 on the right hand side. (The 
remaining terms, e.g. such as the s a on the left hand 
side can be balanced by considering sub-leading correc- 
tions to n(s)). This way of balancing immediately gives 
a = rf/4 and 



A b = v^(s^r d/4 



(73) 



with constraints that d < 2 (Eq. 68 was used to get 
Eq. 71). 

Matching the constant term no on the left hand side of 
(71) is rather counter-intuitive since in the framework of 
Laplace transform constant can normally be disregarded 
when large t behavior is sought for. To see how this 
comes about it is useful to turn back to Eq. (60). 

Equation (60) comes from Eq. (57). For simplicity rea- 
sons we focus on the case A = oo in (57) . It is clear that 
at the right hand side of equation (57) the sn(s) term 
is sub- leading to n . (True enough, n is constant but 
it is multiplied by g(s).) Thus U2(s) indeed has to be 
matched with g(s)n . This procedure results in ampli- 
tude Ab obtained previously. Also, by using form (57), 
one can show that amplitude Ab as given in Eq. (73) is 
valid even for d > 2. Analysis can be repeated with finite 
value of A with the same outcome. Before proceeding, it 
is important to mentioned that procedure outlined above 
does not work at d > 4 and it has to be modified. 

The main findings of this subsection are twofold. First, 
it has been shown that dressed tree calculation is equiv- 
alent to WBGA-I. Second, it was shown that WBGA-I 
(and dressed tree calculation) fail to describe A+A reac- 
tion predicting d/4 decay exponent; 



n(t) ~ V™o( 87r£) *r d/4 



(74) 



All findings of this subsection are summarized in Fig. 1. 
The numerical treatment of Eqs. (54) and (55) confirms 
asymptotic decay given in Eq. (74). Equations (54) and 
(55) were solved previously numerically in ref. [6]. In 
here, the features of the decay curves are somewhat dif- 
ferent from the ones obtained in ref. [6]. For example, 
curves shown in this work have concave form (bend up- 
ward), while curves in fig. 1 of ref. [6] are convex (bend 
downward) as if asymptotics have not yet been reached. 
Also, in here, there is no intersection of curves, which can 
be found in ref. [6] . These difference could come from the 
numerical treatment. Apart from using double precision, 
to obtain curves in Fig. 1 the method of integration was 
used which integrates exactly J* dt'(t — t')~ d / 2 f(t') pro- 
vided f(t) is piecewise linear in t. The more detailed 
description of numerical treatment is shown in the ap- 
pendix B 

In the next subsection it will be shown that, in the case 



of A+A reaction, weaknesses of WBGA-I method extend 
to WBGA-II level. 



B. The WBGA-II approximation 

When term nTk is kept in Eq. (48), equivalent of 
Eq. (55) reads 

<p(t) = n{t) 2 - A / dt'I(t,t')(f(t') (75) 
Jo 

while Eq. (54) stays the same. The I(t, t') is given by 



I(t,t') = G(t-t') exp 



-4A J dt"n(t") 



(76) 



The asymptotics of Eq. (75) can not be extracted by 
Laplace transform, and it is more convenient to use ap- 
proach of ref. [6]. For large t, Eq. (75) can be approxi- 
mated by 



ip(t) *n(() 2 -J((,0)I(i) 



(77) 



where X(t) = A J Q dt'ip(t'). This step is valid provided 
two conditions are satisfied. First, the term I(t, t') has 
to vanish as time difference t — t' grows. Second, the 
integral J Q °° dtip(t) has to be finite. Using (54) one gets 
X(t) = uq — n(t) w no and (77) becomes 

^^-\n 2 + \I(t,Q)no (78) 

Equation above is solved with the assumption that 
asymptotically n(t) <~ A/t, which is checked self consis- 
tently at the end. Using postulated asymptotics for t, one 
can see from Eq. (76) that l(t,0) ~ const t~W 2 + 4XA \ 
Assuming that 



n(t) 2 



, t 



oo 



(79) 



one can solve (78) in form dn/dt = —An 2 , and get 
A = 1/A. Assumption (79) is correct provided 2 < 
d/2 + AXA = d/2 + 4, which is true for any d. This 
shows that n(t) w l/(At) is asymptotic form for the solu- 
tion of Eqs. (54) and (75). This means that the last term 
(nr^) in (48) only influences intermediate behavior when 
t is not too large. For large t, WBGA-II gives exactly the 
same asymptotics as the pure mean field treatment. 

Thus main finding so far is that both WBGA-I and 
WBGA-II fail to describe A+A reaction. This is some- 
what surprising as even simplest pair-approach, e.g. 
Smoluchowskii method, describes exponent of A+A cor- 
rectly. Clearly A+A reaction can not be viewed as a 
weakly interacting Bose-gas. The question is what is the 
minimum modification of WBGA which will provide cor- 
rect result for the A+A model? This question will be 
answered in the next subsection. 
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C. The hybrid of WBGA and Kirkwood 
superposition approximation 

To see how to improve WBGA one has to clarify what 
went wrong in the first place. We start from problematic 
equation (40) which becomes (48) when terms with a£ B 
are set to zero. To trace why WBGA fails it is useful 
to rewrite Eq. (40) as it looks one step before WBGA is 
made, and we keep only terms describing A+A reaction: 



-T k = -2Dk 2 T k - 



a k n 



q=£k 



-f (3) (80) 



where 

f (3) = — \= V a q ((a- k a k - q a q ) + (a k a- k - q a q )) (81) 

is focus of present subsection. In technical terms, the 
usage of WBGA can be translated into approximating 
three point density (a x a y a z ) in a particular way. In the 
case of WBGA-I one simply takes 

(a kl a k2 a k;i ) = (a x a y a z ) = (82) 

while in WBGA-II one assumes 

(a kl a k2 a k3 ) rts 6 kl ,o6 k2 ,o6 k3 ,oa + 

+6 kl ,o6 k2 flS k3fi ao(a k2 a k3 } + 
+6h 1 ,o6h 2 ,o6k 3 ,oao{a>k 1 a>k 3 ) + 
+6k 1 ,o6k 2 ,o6k 3 ,oao(ak 1 ak 2 ) (83) 

where 6k,o = 1 — S k ^- Inserting (83) into (80) gives the 
terms describing A+A process in (40). It is useful to 
transform approximation above into the x-space to un- 
derstand nature of approximation better. Inverse Fourier 
transform of (83) gives [25] 

(a x a y a z ) w n 3 + n((a x a y ) - n 2 ) + 

+n((a x a z ) - n 2 ) + n((a y a z ) - n 2 ) (84) 

The WBGA approximations given in Eqs. (82) and (84) 
can be contrasted to the Kirkwood superposition approx- 
imation, 



(a x a y a z ) w -^(a x a y ) (a x a z ) (a y a z ) 



(85) 



which is known to describe A+A problem correctly (at 
least the decay exponent). 

It is interesting to note that instead of (85) one could 
carry out Kirkwood superposition approximation in a dif- 
ferent way, e.g. as (n x n y n z ) w (n x n y )(n x n z }(n y n z }/n 3 . 
(There is no such ambiguity when only single occupancy 



of site is allowed.) When multiple site occupancy is al- 
lowed, as considered here, we argue that Eq. (85) is more 
accurate way to implement Kirkwood superposition ap- 
proximation. 

By looking at equation (84) it is possible to under- 
stand WBGA better. Eq. (84) suggests that WBGA is 
somewhat equivalent to the additive expansion of corre- 
lation functions. It has been argued that such additive 
approximation is inferior to the Kirkwood superposition 
approximation [4,5]. This in turn explains why WBGA 
fails to describe the A+A reaction. 

Clearly, to correctly describe A+A reaction one has 
to use Kirkwood superposition approximation, which we 
modify further in the spirit of WBGA, assuming that 
higher order products of annihilation operators a k , b k 
with k =/= arc small, and accounting for thermodynamic 
limit (extracting ao or b out of field theoretic averages) . 
This is done in two stages. First, Kirkwood superposi- 
tion approximation in Eq. (85) is rephrased in the k-space 
which gives [26] 



(a kl a k2 a k3 ) 



<5(fci + k 2 + fcjQ 

n 3y3/2 



k 2 +l 



(86) 



To get the improved form for three body terms one in- 
serts (86) into (81) which gives 



f(3) 



n 3 V 2 



(87) 



The expression above was obtained by using symmetry 
properties of T k = T- k and a q = o- q . Also upon insert- 
ing (86) into (81) the two terms on the right hand side 
of (81) contribute equally resulting in factor two in (87). 
Equation above is too complicated to be used in practise 
and we approximate it further. 

In the spirit of WBGA the terms in (86) which con- 
tain large number of correlation functions with fc-vector 
different from zero are neglected. This is done in two 
stages, first sum over q is split into q = k and q =/= k 
parts, and then for each of sums various contributions 
from sum over I arc distilled to extract non-fluctuating 
k = operators. This gives 



f(3) 



i 3 V 2 



[a k (T 2 T k +T r 2 )+ 



+ E cr 9 r (r fc r 9 + r fe r fc _ 9 + r 9 r g _ fe ) 

q^k 



(88) 



where terms of the type T^r^r; with k,q,l 7^ have 
been neglected. Eq. (88) is still complicated to be used, 
at least to obtain analytic result. First term on the right 
hand side of (88) can be neglected since it leads to mean 
field behavior (as shown in section IV B). Second term 
can be absorbed into the one of the three terms under the 
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sum sign, .e.g. under the first term. Second and third 
terms under the sum sign couple correlation functions in 
a non-trivial way and are neglected in the following for 
simplicity reasons. With T w n 2 V, and applying recipe 
just described gives 



(89) 



which is midway between WBGA-II (Eq. 83) and Kirk- 
wood superposition approximation (Eq. 86). It is inter- 
esting to note that equation above is very close to the 
shortened Kirkwood superposition approximation dis- 
cussed in refs. [4,5]. Using approximation above to de- 
couple three body density, and particular form for o^y 
used throughout this section, gives following equations of 
motion, 



^-T k = -2Dk 2 T k - X(n 2 + $) - 2XT k 



(90) 



which should be contrasted with Eq. (48). The most con- 
venient way to solve (47) and (90) is to introduce Xk as 
Tjt = n 2 Xk and n 2 + $ = n 2 x where 



X = l 



(91) 



Xk with k — is set equal to V and does not change in 
time (thermodynamic limit). Applying change of vari- 
ables just described modifies Eq. (47) into 



J^n(t) = -n(t)n{t) 2 

where effective reaction rate 

«(*) - A X (t) 



(92) 



(93) 



was introduced. Same change of variables transforms 
Eq. (90) into 



d 

-Q t Xk = -2Dk 2 Xk - Ax 



(94) 



Also, one can integrate Eq. (92) which gives 

n(t)= U0 , «(*)= fdt'n(t') (97) 
1 + n K{t) J 



It is not possible to obtain closed expression for n{t) and 
n(t). However, asymptotic form of n(t) can be extracted. 

Inserting small s expansion of g(s) (see Eqs. 66 and 
67) into (96) gives 



r(i-d/2) s " < z 



k(s) <~ < 



8ttD 
s[7B+ln(»? s )] 

A „ — 1 

l + Aff(0) S 



d = 2 
d > 2 



(98) 



for s — ► 0. Taking inverse Laplace transform of equation 
above gives leading order behavior for effective reaction 
rate constant, 



«(*) 



r(i-d/2)r(d/2) t " < 2 

871-15 



lnt/rj 
A 

1+Ag(0) 



d = 2 
d > 2 



(99) 



when t — > oo. To find inverse Laplace transform of k(s) 
for d = 2 (second line Eq. 98) is somewhat involved, 
please see appendix C for details. Finally, inserting (99) 
into (97) gives following asymptotics 



n(t) 



r(l-d/2)F(l + d/2)(87rL>i)-< i / 2 d<2 

In SnDt ^ 2 

d > 2 
(100) 



STrDt 

I [i + sW]*- 1 



where g(0) entering in the third row can easily be found 
from (66). 



V. WBGA APPLIED TO THE A+B REACTION 



Equation (94) can be solved for all Xk and k ^ (pre- 
tending that x(t) ls known), and after summing over 
k 7^ one gets following integral equation 



X(t) = l-A f dt'G{t-t')x{t') 
Jo 



(95) 



Solution of the equation above can be found by Laplace 
transform which gives 



k(s) 



X 



s[l + Xg(s)} 



(96) 



Equations of motion for density and correlation func- 
tions describing A+B model result from (35)-(36) and 
(40)-(42) by using (46) with A = 0. For simplicity rea- 
sons, we focus on n^^ng = n case and omit labels A and 
B. Also, as in the previous section n = n ,A = no.B- 
Applying procedure outlined above leads to the equations 
for particle densities 



where $ c is given by the AB correlation function, 



(101) 
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(102) 



Equations for correlators T k = T AA = T B B and = 
T AB are given by 



- + 2Dk 2 ) T c k = -S(n 2 + $ c ) - 2Sn(T k + T c k ) (103) 



§- t + 2Dk 2 ) r,, = - 2rt//(r A . + r' A ,) 



(104) 



The most convenient way to solve equations above is to 
diagonalise them by subtraction and addition. The final 
result is that correlations of AB pairs are governed by 

#c(*) = S f dt'[G(t, t') + I(t, t')][n(t') 2 + $ c (0] 
Jo 

(105) 

and for AA (or BB) pairs as 

$(i) = <S /*dt'[G(t,t') -■?(*, OIK* 7 ) 2 + *c(*')] ( 106 ) 
Jo 

The J(t, i') appearing in equations (105) and (106) has 
the same form as in (76) with trivial change of A into 5. 
Same type of analysis as in the subsection IV B leads to 
conclusion that approximation G(t, t') + I(t, t') w G(t, t') 
can be used, which in turn leads to d/A density de- 
cay exponent. Interestingly enough, both WBGA-I and 
WBGA-II approaches lead to correct d/4 exponent when 
used to solve A+B model. 

In some sense the WBGA approach seems to be suited 
rather well for A+B reaction. Quite the contrary can be 
said for A+A reaction as shown previously. To under- 
stand workings of WBGA on the more general model the 
ABBA model will be studied in the next section. 

VI. WBGA APPLIED TO THE ABBA MODEL 



The ABBA model was suggested in refs. [14,15] and 
it has a number of interesting properties. The model is 
obtained as the special case of the general two species 
model discussed in section II where A+A and B+B reac- 
tions occur with same reaction rate A, while A+B with 
reactions rate 5, and S > A is taken. The fact that A and 
B species have equal diffusion constants is very impor- 
tant since it leads to survival of minority species. Minor- 
ity species is the one with smaller concentration at t = 0, 
and we chose B to be minority, i.e. n ,A > no,B- By 
survival it is meant that particle density ratio saturates 
to constant u(t) = nA(t)/riB{t) — ► const as t — > oo. (The 
mean field analysis predicts vanishing of minority species 



u{t) — > oo as t — > oo. Also, similar model has been stud- 
ied in ref. [27] where it was shown that if Da ^ D B only 
one type of species survives.) 

The goal of the present section is to understand 
strengths and weaknesses of WBGA approach by testing 
it on the ABBA model which combines A+A and A+B 
reactions. These reactions have been studied in previous 
sections but now they are allowed to occur simultane- 
ously. It is interesting to see how WBGA approach per- 
forms in such situation. Clearly, there are plenty of possi- 
bilities how to combine A+A and A+B reactions, but the 
way they are combined in the ABBA model gives rise to 
many interesting effects. (For example, there is survival 
of minority species, recovery of initial density ratio, de- 
pendence of amplitudes on initial densities and reaction 
rates [28], please see refs. [14,15] for details). 

We start from equations of motion given in section III 
which describe very general two species reaction-diffusion 
model. Assumptions in Eq. (46) describe content of the 
ABBA model. Using (46) in Eqs. (35)-(36) gives 



—n A = - (Xn 2 A + Sn A n B + \$aa + 5$ab) 



dl nB 



where 



(Xn 2 B + 6n A n B + A$ BB + 6<$>ab) 



p,u 



A,B 



(107) 
(108) 

(109) 



Also, Eqs. for correlators (40)-(42) simplify to 

= -2Dk 2 T£ A - A (n A + *aa) 
-2(2\n A + Sn B )T AA - 25n A T AB 

d 



gT BB = -2Dk 2 Y BB A (n% + $ SB ) - 

-2(2An B + Sn A )T BB - 2Sn B T A 
^-T AB = -2Dk 2 T AB - 5 (n A n B + <& AB ) 



-.AB 



-(2\ + 6)(n A + n B )T AB - 
-5(n A T BB + n B T AA 



(110) 
(111) 

(112) 



The equations above will be solved in the next two sub- 
sections within WBGA-I and WBGA-II approaches. 



A. The WBGA-I approximation 



In the framework of WBGA-I all seemingly 0(n 3 ) 
terms of the type n p T^, with p = A, B and v = 
AA,AB,BB are thrown away in (110)-(112). Following 
similar steps as in subsection IV A gives 
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gj. n A = ~{\<PAA + StfiAB) 
^"B = ~(^<PBB + 5lf AB ) 



and 



y AA =n 2 A -\ f dt'G(t - t')cp AA {t') 
Jo 

<PBB=n 2 B -\ [ dt'G{t - t')ip BB (t') 



VAB = n A n B 



5 / dt'G{t-t')y AB {t') 



(113) 
(114) 

(115) 
(116) 
(117) 



where (p pl/ = n p n v + <f> p „ for p, v e {A, B}. To solve 
Eqs. (113)-(116) it is possible to employs same technique 
as in section IV B. Equations above can be approximated 

by 



0^n A -G{t 7 0)l AA (t) 
0^n 2 B -G(t,0)l BB (t) 
w n A n B - G(t,0)l AB (t) 



(118) 
(119) 
(120) 



where 



T pv {t) = (M p , v + 55 p , v ) f dt'<p pv {t') , P = A,B (121) 
Jo 

By integrating Eqs. (113) and (114) a useful relationship 
can be derived for 2" p „, p, v = A, B; 

T-aa (t) +I AB (t) = n , a - n A {t) w n ,A (122) 
lBB{t) +l AB {t) = n ,B - n B (t) w n ,s (123) 



Using (118)-(120) gives (n A +n B ) 
2T AB ), and adding (122) and (123) finally gives 

(n A + n B f w G(t)(n ,A + no,s) 



G(t)(X AA +X BB + 
(124) 



Also from (118) and (119) n\ - n| w G(t)(J AA - X BB ) 
and subtraction of (122) and (123) gives 



- n| w G(t)(n ,A - n ,s) 



(125) 



After solving (124) and (125) for n A and n B one gets 



V n o : A + riQ.B 



(8irDty 



A,B 



(126) 



According to WBGA-I both particles decay with d/A 
exponent and amplitudes given above. The WBGA-I pre- 
dicts same decay exponent as for the pure A+A model. 
As in the case of A+A reaction, the value for d/A expo- 
nent obtained here is not correct. The computer simu- 
lation and e-expansion analysis of this reaction suggest 
d/2 exponent [14,15]. To see what happens when 0(n 3 ) 
terms are kept in (110)-(112) we proceed with WBGA-II 
calculation. 



B. The WBGA-II approximation 



In the WBGA-II approximation, when all terms are 
kept in Eqs. (110)-(112), it is useful to rewrite these equa- 
tions in the vector form 



d_ 

di 



2Dt 











( XipAA \ 


(IP 


= -P(t) 
















V 8<PAB J 



(127) 



where the matrix P is given by 

4A + 2Su 25 
P = n A | 4Au + 25 25u \ (128) 

5u 5 (2A + 5)(l + u) 

with u — n B /n A . 

Vector equation (127) is very hard to solve analytically. 
However, there are some guidelines how to extract late 
time asymptotics. At the WBGA-I level it appears that 
ABBA model and the A+A model are very similar. In 
the following it will be assumed that such similarity can 
be extrapolated to the presently studied WBGA-II level. 
This implies that mean field behavior should be expected 
from Eq. (127). 

To get the feeling for what follows it is useful to analyze 
Eqs. (107) and (108) at the mean field level by neglect- 
ing fluctuations and setting <fr pl/ to zero with p, v = A, B. 
Carrying out such procedure gives 



d_ 

di 

d_ 

di 



n A 



n B 



(\n 2 A + 5n A n B ) 
(An| + 5n A n B ) 



(129) 
(130) 



Equations above can be solved approximatively for large 
t (please see refs. [15] for details) and one obtains 



n A 
n B 



1 

At 



no,B 



[no, A Xt}" 



(131) 
(132) 



provided 7 = 5/X > 1 and no, A > no lB . For 5 = A 
(7 = 1) or n A = n , B solution is trivial, and it can be 
easily shown that in such case the ABBA model belongs 
to the A+A universality class. These simple cases are not 
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considered here. Initial imbalance in particle concentra- 
tion leads to faster diminishing of minority species, i.e. 
u(t) — nB{t)/riA{t) — > as t — ► oo given < u(0) < 1. 

In the following we assume the mean field ansatz (131)- 
(132) and try to solve (127) with it. The validity of 
such mean field ansatz will be checked self consistently at 
the end. For large times, and with mean field behavior 
(u — ► 0), the matrix P can be approximated by 



with 



n 



4 2 7 
2 7 

7 2 + 7 



(133) 



(134) 



The fact that P (in the approximate form) is constant 
matrix multiplied by time dependent function implies 
that [P(i),P(t')] = (dot over symbol P denotes time 
derivative). This being the case, Eq. (127) can be treated 
as scalar equation and calculation similar to the one in 
the subsection IV B gives 



/ <PAA{t) \ 

<PBB(t) 
V <PAB(t) J 



n A (t) 2 
n B (t) 2 
\ n A (t)n B (t) ) 



f 

Jo 



dt'3(t,t') 



Xlfi BB (t') 



(135) 



(Ji = 4 , u} 2 = 27 , lu 3 = 2 + 7 



and matrix U contains eigenvectors 



(139) 




(140) 



Inserting (131) into (137), and assuming large t, leads to 
f(i,0) ~ const + In t (141) 

and using (141) in (136) gives 







J(t,0) -const t~ d/2 V [ t- 2 '< j U- 1 

V H 2 +t) J 

(142) 

Finally, the second term on the right hand side of (138) 
can be calculated explicitly. Inserting (142) into (138), 
and assuming that I pv p,v = A, B are constants (can be 
checked for self consistency at the end), results in 



<fiAA ^n\+ t- d ' 2 { Cl t- Ul + c 2 i-<" 2 + cat'" 3 ) 

fBB ~n 2 B + t- d ' 2 c i t-° J2 

VAB w n A n B + t- d,2 (c 5 t- U2 + c 6 i _W3 ) 



(143) 
(144) 
(145) 



where matrix 3(t,t') is given by 

3{t, t') = G(t, t')exp [-£(t, t')U] (136) 

and 

£(t,t') = A / dt"n A (t") (137) 
Jv 

Please compare Eqs. (75)-(76) and (135)-(136). They 
are very similar, the only difference being in the matrix 
character of (135-136). Following same steps as in section 
IV B the Eq. (135) can be approximated as 



/ VAA ^ 




( < \ 




(Iaa \ 


<Pbb 




n% 


-J(t,0) 


Tbb 


V <Pab J 




V n A n B ) 




\Tab) 



Now we proceed to show that, as in the case of (77), the 
second term on the right hand side of equation (138) can 
be neglected. 

Matrix II can be diagonalized as IIU = UJ7. The ft 

is diagonal matrix containing eigenvalues 



The explicit form of constants C\, C2, c 3 , C4, C5 and Cq 
is not interesting since aim is to show that terms con- 
taining these constants are sub-leading to the mean field 
terms. By studying equation above row by row, it is pos- 
sible to show that for 7 > 1 terms involving constants 
are sub-leading to the mean field terms. 

To see that terms originating from 3(t, 0) are sub- 
leading one really has to calculate U explicitly. For 
example, not knowing that contribution from u>i is ab- 
sent in (144), there would be a need to compare t~ 21 
(asymptotics of the mean field n 2 B term in 144) with 
£-(d/2+4) (coming from 3(t, 0) and ui\ eigenvalue). One 
would conclude that 7 = 5/X can not be too large if 
mean field asymptotics is to hold. In reality, there is no 
such bound on ratio 5/X since eigenvalue u\ does not ap- 
pear in Eq. (144), but this can only be seen after explicit 
calculation. 

Main findings so far is that WBGA describes ABBA 
and A+A model in the same way. For both models the 
WBGA-I (WBGA-II) predict d/4 (mean field) density 
decay exponents. In the next section attempt will be 
made to improve WBGA method in order to obtain cor- 
rect value of density decay exponent for ABBA model. 
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C. The hybrid of WBGA/Kirkwood applied to A+A 
and B+B sectors 



with p,v = A, B. Using (150) and (151) in Eqs. (113) 
and (114) results in 



How deep the weakness of WBGA goes? What needs 
to be changed in equations of motion (110)-(112) in order 
to get correct decay exponent? To answer these ques- 
tions we begin by modifying the A+A and B+B reaction 
sectors in (110)-(112) by using the recipe from section 
IV C. It was already remarked in the section II that 
contributions to H describing different reaction sectors 
enter additively, and this feature is reflected in equations 
of motion (110)-(112). Because of this it is possible to 
focus on terms describing influence of A+A and B+B re- 
actions, i.e. terms proportional to A, which govern decay 
exponents of the ABBA model. At the moment terms 
proportional to 8 in (110)-(112) that describe A+B reac- 
tion will not be changed. 

If one is to follow procedure described in section IV C, 
the 0(n 3 ) term in Eq. (110) has to be modified as 



4\n A T AA - 2Ar£ 



and likewise for Eq. (Ill) 



AA n A + ^AA 



n A 



(146) 



d 

-Q t n A = -(^n 2 A XAA + SuaubXab) (152) 
d 

-Q-riB = -(An|xss + Sn A n B XAB) 



(153) 



Implementing same notation in Eqs. (148), (149) and 
(112) gives 



d_ 

dt 

d_ 

dt 

d_ 

dt 



xt A = 



x BB = 



x AB = 



-2Dk 2 X AA ~ ^XAA - 

-25n B (l- X AB)x AA ~2Sn BX £ B 

-2Dk 2 X BB - Xxbb - 

-2Sn A (l - XAB ) X BB -26n AX AB 

-2Dk\ AB - 8 X ab - 

-[Xn A (2 - xaa) + An B (2 - xbb) + 

+5(n A + n B )(l-XA B )}xi B - 
-5(n AX £ A + n BX BB ) 



(154) 



(155) 



(156) 



A\n B T BB 



2Ar£ 



BB il b 



BB 



n B 



(147) 



This gives a new set of, hopefully better, equations: 

a. 



r AA 



2Dk 2 T AA - A (n 2 A + $ AA ) 



^AA 11 A 



AA 



n A 

AA i „ r AB\ 



dt 



r BB _ 
1 fc — — 



-2\T k 

-2S(n B T£ A + n A T A 
2Dk 2 Yl B - A (n% + $> BB ) - 
BB n 2 B + <5>BB 



(148) 



-2Arf 



n B 

^BB i „ rAB\ 



(149) 



Equation (112) stays the same, although the Eq. (112) 
contains term proportional to A which should be modi- 
fied if one follows the principle outlined above. However, 
at the moment, Eq. (112) will not be changed. 

It is useful to employ similar notation to the one used 
in section IV B; 



rnv pv 



and 



Xpv 



\_ \ ^ pu 
1/ 2-^1 *k 



v 



(150) 



(151) 



fc#0 



The numerical solution of the set of equations above 
is shown in Fig. 2 (dotted line). The full line is a re- 
sult of Monte Carlo simulation where particle densities 
are obtained as ensemble averages over 500 runs (simu- 
lation is repeated 500 times with a shift in the random 
number generator). The Eqs. (154)-(156) do not describe 
ABBA model properly, not even qualitatively, since mi- 
nority species die out faster (the particle density ratio 
grows to infinity) while the simulation shows that den- 
sity ratio should saturate to a constant value (full line). 
Thus the attempt of modifying terms describing A+A 
and B+B reaction sector in AA and BB correlation func- 
tions by using shortened Kirkwood superposition approx- 
imation does not cure weaknesses of WBGA approach. 
The inspection of individual density decays reveals that 
the equations above correctly describe decay of majority 
species, i.e. n A <~ (const)i~ d / 2 , but fail do describe decay 
of minority species n B . 

In the following more terms will be modified by using 
Kirkwood superposition approximation. The dash-dot 
line in Fig. 2 shows a solution of equations (not shown 
here) obtained from modifying A+A and B+B reaction 
sectors using shortened Kirkwood superposition approx- 
imation in the Eq. (112) describing time evolution of 
AB correlation function. Equations obtained in this way 
are identical to the (154)-(156) with only difference that 
Eq. (156) changes in a way that terms proportional to A 
drop out. In Fig. 2 it can be seen that even if all A propor- 
tional terms are modified the density ratio curve climbs 
to infinity, which indicates that minority species dies out, 
which is not correct behavior. However, the overall trend 



15 



gets better as the dash-dot curve lies bellow dotted one 
and is pushed towards simulation curve. 

To continue this line of incremental changes the equa- 
tions of motion where studied where even the A+B re- 
action sector was modified using shortened Kirkwood 
superposition approximation (all terms proportional to 
5 were modified). Equations obtained in this way are 
same as in (154)-(156) the only difference being in the 
fact that all seemingly 0(n) terms drop out. Thus in 
the (154)-(156) only diffusion term and terms X\aa, 
^Xbb and Sxab are kept in (154), (155) and (156) re- 
spectively. These equations are not shown explicitly to 
save the space but it should be clear how they look 
like. The numerical solution for this set of equations 
is shown in Fig. 2 as dashed line. [It is possible to ana- 
lytically extract density decay asymptotics for this trun- 
cated set of equations which gives n^(t) ~ (const)t _d / 2 
and n B (t) ~ (const')^/ 4 .] 

The set of equations where Kirkwood superposition ap- 
proximation has been applied fully agrees with the nu- 
merical simulation much better than the rest which are 
mixture of WBGA and Kirkwood superposition approx- 
imation. This is strong indication that, at least for the 
ABBA model, the Kirkwood superposition approxima- 
tion is superior to the WBGA method. For example, in 
Fig. 2, the trend in all curves improve as the content of 
superposition approximation is increased in equations of 
motion (dotted line is worse as it climbs fastest, dash-dot 
line is a bit better, while only dashed line where super- 
position approximation is implemented fully saturates to 
constant). As the goal of the present study is to under- 
stand WBGA method better, the more detailed analysis 
of ABBA reaction based on Kirkwood superposition ap- 
proximation will be presented in forthcoming publication. 



VII. CONCLUSIONS 

The goal of this study was to test workings of WBGA 
and to relate some seemingly independent calculation 
schemes available in the literature. In particular, this 
work has addressed few important issues. 

(1) It was shown that WBGA-I is equivalent to the 
dressed tree calculation introduced in ref. [19], and this 
equivalence holds for any model where particles annihi- 
late in pairs. Furthermore, it was shown that for the 
A+A reaction the dressed tree calculation results in d/A 
density decay exponent. Thus the contradictory claims 
of refs. [19] and [6] have been sorted out. 

(2) The WBGA method describes A+B reaction well, 
but does not work for A+A and ABBA models, and 
it is reasonable to expect that there are more models 
that can not be described by WBGA. In the case of 
ABBA model WBGA predicts faster vanishing of mi- 
nority species which is suggestive of the A+B like be- 



havior rather than the behavior of the ABBA model 
as found in refs. [14,15]. This bias towards A+B type 
behavior is very hard to get rid off as successively cor- 
recting more and more terms in equations of motion for 
ABBA model by using Kirkwood superposition approx- 
imation results in faster vanishing of minority species. 
The vanishing of minority species persists until all terms 
are modified by Kirkwood superposition approximation. 
The WBGA emphasize A+B reaction sector too strongly 
in the ABBA model. 

(3) Findings of this work suggest that formalism em- 
ployed by Mattis-Glasser in [6] where small no expan- 
sion is introduced (and applied to study A+B model) is 
somewhat questionable. This procedure works on A+B 
model, but might not work for other models. It can be 
shown (by rescaling a^no — > a 1 and a/no - > a ) that for 
the type of models studied in here, the small no approx- 
imation of ref. [6] amounts to taking away three body 
terms in Hamiltonian given in Eq. (21) or (28) (e.g. op- 
erators of the type a^aa and likewise any mixture of a or b 
operators). The neglect of these terms amounts exactly 
to WBGA-I approach, i.e. neglect of seemingly C(n 3 ) 
terms in equations of motion for correlation functions; 
for example, equations (54)-(55) of section IV A (A+A), 
or equations (113)-(117) from section VIA (ABBA). In 
here it has been shown that WBGA-I set of equations re- 
sult in incorrect d/A density decay exponent when used 
for A+A and ABBA models. For these reasons, small 
no expansion, which effectively means taking away three 
body term in Hamiltonian, can not be trusted if used 
beyond A+B model. 

(4) The Kirkwood superposition approximation was 
formulated for the case of lattice dynamics with multiple 
occupancy of lattice sites allowed. There is strong indi- 
cation from the present analysis that Kirkwood superpo- 
sition approximation is superior to the WBGA method, 
at least when applied to the A+A and ABBA models. 
However, one has to be careful in claiming superiority 
of Kirkwood superposition approximation over WBGA 
approach since in here it was shown that Kirkwood su- 
perposition approximation describes ABBA model quali- 
tatively, but it remains to be seen whether it works on the 
A+B model which WBGA describes well. Similar study 
(where only single occupancy of lattice sites was allowed 
and reaction range was assumed short but finite) showed 
that Kirkwood superposition approximation can describe 
A+B reaction [4] , and claim of superposition approxima- 
tion superiority is likely to be correct but, nevertheless, 
such claim has to be tested throughly. 

The application of Kirkwood superposition approxi- 
mation to the ABBA and A+B models will be presented 
in the forthcoming publication as there are many inter- 
esting technical details that need to be sorted out. For 
example, from present study there are some indications 
(results not shown here) that when using Kirkwood su- 
perposition approximation for extremely short range re- 
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action, such as on-site reaction studied here, care has to 
be taken so that thermodynamic limit is accounted for 
properly. It seems that one has to approach Kirkwood su- 
perposition approximation through formalism developed 
in section IV C. 

To conclude, it would be interesting to have a rela- 
tively simple approximation at hand, not far away from 
pair approximation, which could be used to extract quali- 
tative asymptotics for arbitrary reaction-diffusion model. 
Clearly such program is ambitious since in reality one is 
bound to make approximation which is related to the 
particular model but, nevertheless, it is worth a try. The 
A+A and A+B reaction-diffusion models (or combina- 
tion of them) are excellent bench-mark models and any 
successful approximation should strive to describe these 
reactions properly. Present study is a attempt in this 
direction. 



APPENDIX B: INTEGRATION OF SINGULAR 
KERNEL 



In here the numerical treatment of Eqs. (54-55) is de- 
scribed in a more detail. The general procedure for inte- 
grating expressions of the type 

h[f] = I* dsK(t h s)f(s) (Bl) 
Jo 

where K(t, s) is singular when s approaches t is described 
in ref. [29]. The i = 0, 1, 2, . . . and ti — ih. In here the 
particular focus in on the particular form of singular ker- 
nel originating from propagator G(t, t'); 

K(t,s)^{t-s + r 1 )- a (B2) 
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The aim is to find quadrature formula which can inte- 
grate (Bl) exactly if /(s) is linear within intervals U-\ , ti. 
Thus goal is to find coefficients Wij in 



(B3) 



APPENDIX A: APPROXIMATION FOR G(t t') 



In Eq. (52) sum over k can be approximated as integral 



such that above equation turns equality for piecewise lin- 
ear function f(s). Implementing this procedure gives 



7 *[/] = E I''** dsK(ti,8)f(s) 

A—n J S-i 



(B4) 



j=0 " s i 



G(t - 1') = 



(2tt) 



dpe 



-2Dp 2 (t-t') 



(Al) and after approximating f(s) as linear function in inter- 
vals tj,tj+i for j = 0, .., i — 1 



The equality sign is strictly valid in thermodynamic limit 
only when V — > oo has been taken. It is useful to ap- 
proximate integral above further by changing bounds of 
integration 



G(t - t') 



1 



(27t)° 



J — ( 



dpe 



-p 2 A- 2 -2DK, 2 (t-t') 



(A2) 



A corresponds to a large cutoff for the fc-vector sum- 
mation (it regularizes UV divergences of theory). The 
comparison of (Al) and (A2) quickly reveal that A is a 
constant roughly equal to ir. The evaluation of integral 
above gives Eq. (53) with rj = 2 p A i ■ It is convenient to 
chose 77 = g^jj and A = 2y / 7r since this gives correct value 
for G(0) = 1. For large t — t' particular value of A ap- 
pears to be irrelevant, i.e. absent from final expressions 
for particle density for large times. However, this really 
depends on the critical dimension of the field theory. 



/(*) 



gives (B3) with 



w i0 = / dsK(U, s) — 
Jo " 



Wa = 



r dS K(t l}S ) s 



h 



+ 1 dsK(t h s) Sj+1 3 



/Sj + l 
j 

f dsK{ti,s) 



h 

S - Sj-i 

h 



(B5) 



(B6) 

(B7) 
(B8) 



and after performing integrals above with K(t, s) given 
in (B2) 
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w l0 



(a-l)(a - 2) 
+ {i + i 1 -l) 2 - a 
h x - a 



[(2 - a - 1 - n)(i + iiy- a + 



[(i-j+r,-iy- a + 



11 - (a-l)(a-2) 1 

+ (i-3+V + l) 2 - a -2(i-3+vf- a } 



(a-l)(a-2) 
+r ] 1 - a (a-rj-2 



[(! + »?) 



2-a 



+ 



(B9) 



(BIO) 



(Bll) 



The pair of equations in (54-55) is discretized as follows. 
First the differential equation (54) is rewritten in inte- 
gral form as nit) = no — A J„ dstp(s) and trapezoidal rule 
is used to evaluate integral since all functions are well 
behaved. However, for Eq. (55) the rule (B3) and (B9)- 
(Bll) designed for singular kernel is used. Implementa- 
tion of this philosophy gives 



i-l 

m = n - Xh[— + 2^f 3 + —\ 

3 = 1 

i-1 

<Pi = nf - w ij<Pj ~ ^ w ii¥>i 

3=0 



(B12) 
(B13) 



where n, — n(ti) and ipi = (p(U) for i = 0, 1, 2, . . .. Given 
that all tij and ipj are known for j = 0, 1, 2, . . . , i — 1 
using equations above it is possible to calculate rij and 
ipi. The iteration is started with n = n(0) and <fo = n o- 



APPENDIX C: FINDING INVERSE LAPLACE 
TRANSFORM OF k(s) FOR d = 2 

In here inverse Laplace transform of k(s) for = 2 
given in Eq. (98) will be found. Due to the presence of 
log one has to use Bramowitz contour to perform inte- 
gration over s. Also, function k(s) does not have poles. 
This means that only contribution to n(t) comes from 
the branch cut and one obtains 



, . n „ f°° du 
nit) = 8ttD / — e 
Jo u 



1 



(■je + ln(r)u)) 2 + n 2 



(CI) 
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which, upon using the fact that for particular range of 
integration above l/u < 1/c and performing remaining 
integration, gives 
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It will be shown that n a (t) vanishes lot slower than Kb {t) . 
One can approximate expression for K a as 
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and it can be shown that terms omitted do not influence 
leading order behavior for K a . By using partial integra- 
tion expression above becomes 
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The integral over u is most conveniently calculated by 
changing variables tu = v which gives 
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By sending upper integration limit to infinity, and ex- 
panding denominator in series over In vj lnt one gets 



8ttL> 
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(C7) 



Keeping in mind that transformation t — > t/n has to 
be made in the expression above gives result for in 
Eq. (99). 



It is useful to re-scale integration variable as rju — > u 
with only change that time variable appears in combina- 
tion t = t/n. In the following we set n = 1 but keep in 
mind that at the end of calculation t has to be changed 
into t/n. 

As t grows, due to the presence of exp(— ut), only 
smaller and smaller values for u contribute to the integral 
above and it is useful to split integration over u in the 
two intervals (a) from to c and (b) from c to infinity, 
where < c < 1; «(t) = K a (t) + Kb(t). It is possible to 
find upper bound for Kb as follows. Starting from 
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FIG. 1. FIG 1. The numerical solution of Eqs. (54-55) 
for d = 1,1.5,2.5,3 (solid lines). The dotted lines indicate 
asymptotics as given by Eq. (74). Time is given in seconds 
and particle density n(t) is dimensionless in units of particles 
per site. Initial density n was set equal to 1, and reaction 
rate A = Is -1 was used. 
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FIG. 2. FIG 2. The numerical solution of Eqs. (152-156) 
for d = 1 with increasing amount of Kirkwood superposition 
approximation embedded (dotted, dash-dot and dashed line) . 
Full curve is result of a Monte Carlo simulation (average of 500 
runs). Parameters used are L = 1000, ua(0) = 2, ns(O) = 1, 
A = 1 and S — 2. 
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